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ABSTRACT 


This  report  explores  methods  for  obtaining  transient 
solutions  to  queueing  problems  which  can  be  represented  in 
the  form  of  differential-difference  equations.  Six  distinct 
methods,  representing  the  most  frequently-encountered  in  the 
open  literature,  are  discussed  as  to  their  value  in 
numerical  work.  The  method  of  Runge-Kutta  integration  of 
these  equations  was  found  to  be  superior  to  the  numerical 
evaluation  of  analytic  solutions  of  a  particular  queueing 
model. 

A  generalized,  Runge-Kutta  programming  package,  written 
in  FORTRAN  IV  for  the  IBM  360/65,  is  presented  and  described 
in  detail  for  use  on  queueing  problems.  Generality  is 
achieved  by  requiring  the  user  to  write  a  subroutine  to 
evaluate  his  queueing  equations  when  required  by  the 
programming  package.  Several  sample  problems  are  presented 
and  solved  to  demonstrate  the  wide  potential  applicability 
of  this  method  and  associated  computer  program.  The 
advantage  of  this  method  is  that  the  usual  rate  parameters 
to  describe  state  transitions  may  be  any  functions  freely- 
chosen  by  the  user. 

The  appendix  which  contains  the  FORTRAN  program  and 
instructions  for  its  use  may  be  separated  fr<\,i  v*'*:  report 
and  used  separately.  ' 
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CHAPTER  1 
INTRODUCTION 

This  thesis  explores  problems  associated  with 
obtaining  the  transient  solutions  of  queueing  problems  with 
varying  rate  parameters  and  various  solution  methods  that 
have  been  employed.  The  work  of  Bhat  (2)  and  Neuts  (11)  has 
shown  that  the  lack  of  numerical  work  in  the  solution  of 
transient  queueing  models  has  "handicapped  practitioners  to 
obtain  meaningful  results'*  (2).  This  thesis  addresses  that 
problem  by  presenting  a  complete,  well-tested  computer 
'■rogramrning  package  for  obtaining  Runge-Kutta  solutions  to 
transient  queueing  problems.  The  thesis  concludes  that 
Runge-Kutta  integration  of  the  differential-difference 
equations  of  queueing  problems  is  "best"  in  several 
respects,  given  the  current  state-of-the-art  of  available 
computational  facilities. 

Chapter  two  deals  with  a  review  of  the  computational 
aspects  of  several  techniques  contained  in  the  literature. 
These  include  closed-form  analytic  solutions, 
approximations,  simulation  and  Runge-Kutta  integration.  In 
Chapter  three,  an  evaluation  of  three  basic  methods  is  made 
on  the  basis  of  computer  programs  which  were  written,  tested 
and  analyzed  for  each.  This  Chapter  utilizes  the  single¬ 
channel,  single-server  queue  as  the  principle  comparative 


model.  Examples  of  solutions  to  other  queueing  models  are 
then  presented  in  Chapter  four.  The  thesis  is  then 
summarized  in  Chapter  five,  followed  by  Appendicies  A,  B  and 
C  which  present  listings  of  computer  programs  used  in 
Chapters  three  and  four. 

The  programming  package  for  Runge-Kutta  solution  of 
queueing  problems  is  presented  in  detail  in  Appendix  C. 
Included  with  it  are  complete  instructions  to  guide  users  in 
its  application  to  their  problems  and  a  detailed  explanation 
of  a  sample  problem. 
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CHAPTER  2 

REVIEW  OF  METHODS  FOR  PARTICULAR  TRANSIENT  QUEUES 

Queueing  problems  are  frequently  modeled  in  the  form  of 
a  set  of  simultaneous,  first-order,  differential  eauations 
called  the1  Chapman-Kolmorgorov  (5)  equations.  These  are 

P(t)  =  AM t)  , 

where  P(t)  and  P(t)  are  column  vectors  and  A  =  {a. .}  is  a 
-  - 

system  transition  matrix  representing  the  time-rate  of  a 

transition  from  state  i  to  state  j.  If  P(t)  =  [PQ(t), 

P  (t),***,  P^ft)]  ,  then  pn(t)  represents  the  probability 

that  the  system  is  in  state  n  at  time  t.  Correspondingly, 

F»n ( t)  represents  the  time  rate  of  change  of  P^ft).  In 

queueing,  the  system  is  usually  thought  of  as  a  set  of 

channels  and  servers,  related  by  a  queue  discipline  which 

describes  the  flow  of  items  in  the  system.  Thus  the 

notation  P  (t)  describes  the  time  fluctuation  of  the 
n 

probability  of  seeing  n  items  in  the  system  for  0  <_  n  <  °°. 

For  the  M/M/1  queue  (15),  the  general  equation  for 
n  £  1  in  the  system  is 

=  >Pn-1(t)  ‘  (UulPn(t)  +  UPn+1(t)- 
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The  resulting  matrix  form,  limiting  the  maximum 
number  in  the  system  to  N  is 
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In  the  transient  queueing  problems  considered  by  this 
thesis  it  is  desired  to  admit  rate  parameters  a^j  of  any 
arbitrary  form  and  functional  dependence  on  either  time  t  or 
number  in  the  system  n.  Many  analytic  techniques  require 
that  n  tend  to  infinity  for  a  solution.  In  this  thesis,  the 
Runge-Kutta  methods  discussed  will  require  n  be  definitely 
bounded  at  some  finite  value  N.  Further  consideration 
regarding  this  point  will  be  presented  in  Appendix  C. 

All  too  frequently,  only  steady-state  results  are 
considered,  i.e.,  the  solution  of 


0  =  A  P 
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even  when  the  system  under  study  has  frequent  inherent 
changes  in  arrival  or  service  rates.  Numerous  situations  can 
be  found  where  arrival  rates,  service  rates  and  perhaps 
other  system  parameters  are  changing  continually.  One  such 
situation  is  in  a  repair  shop.  Initially,  the  shop  is  empty 
at  the  start  of  the  day.  Arrivals  occur  at  a  rate  greater 
than  the  service  rate  during  the  day.  Steady  state  solutions 
to  such  systems  are  inadequate.  What  is  needed  is  a  method 
which  will  yield  transient  solutions  to  such  problems  with 
no  more  difficulty  than  that  required  for  steady  state 
solutions.  That  method  is  presented  in  this  thesis. 

The  argument  for  seeking  a  standard  technique  for 
solving  a  wide  class  of  transient  queueing  problems  has  a 
strong  practical  appeal.  The  technique  can  be  studied  on 
problems  with  known  solutions  in  order  to  determine  its 
properties  under  wide  ranging  conditions.  The  analyst  can 
thus  develop  intuitive  insights  which  are  helpful  when 
applying  the  technique  to  problems  whose  solutions  are  not 
known.  In  applying  a  standard  solution  technique,  however, 
the  analyst  runs  the  risk  of  overlooking  a  possible  analytic 
solution  had  he  tried  some  other  methods.  In  Chapter  three 
the  economics  of  attempting  to  solve  the  analytic  closed 
form  of  a  transient  queueing  problem  versus  using  Runge- 
Kutta  and  Monte-Carlo  techniques  will  be  discussed  further, 

A  search  of  the  literature  reveals  that  transient 
solutions  to  numerous  queueing  problems  can  be  found. 
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However,  in  most  instances  these  solutions  appear  in 
graphical  form  with  no  reference  to  the  numerical  techniques 
used  to  obtain  them.  In  cases  where  the  numerical  techniques 
employed  are  revealed  (12,  13)  little  discussion  has  dealt 
with  the  reason  for  selecting  the  particular  technique  from 
among  the  many  available.  Therefore,  the  literature 
concerning  explicit  techniques  for  transient  queues  is  quite 
limited  although  at  least  eight  techniques  have  appeared  in 
the  open  literature.  Six  of  the  most  useful  are 

1.  Use  of  generating  functions  (z-transforms) 

2.  Analytic  continuation  of  the  differential  equations 

3.  Series  expansion  of  rate  parameters 

4.  Transformation  to  integral  equations 

5.  flonte-Carlo  simulation 

6.  Direct  integration  of  the  differential  equations 

Two  additional  methods  do  not  appear  promising  for 

study  here.  These  are  the  method  of  numerical  inversion  of 
the  Laplace  transform  and  the  method  of  eigenvectors.  The 
method  of  numerical  inversion  does  not  offer  the  desired 
qualities  of  generality  and  numerical  stabilty.  The  Laplace 
transform  is  difficult  to  obtain  when  rate  parameters  are 
functions  of  time.  Also,  in  a  thesis  by  Schlenker  (16, 
p. 1 4)  the  numerical  method  proposed  by  Bellman  is  shown  to 
be  frequently  unstable  and  limited  in  accuracy. 


The  method 

of  eigenvectors 

also 

does 

not 

promising.  In  a 

paper  by  Kaspar 

(8)  , 

the 

method 
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discussed  on  a  multiple,  parallel  channel  system  with  no 
queues  allowed.  Kaspar  claims  that  if  the  A  matrix  can  be 
shown  to  have  distinct  eiqenvalues,  then  a  new  system  of 
linear  differential  equations  with  separated  variables  can 
be  obtained.  In  our  problem,  however,  the  matrix  A  may  be 
time-varying  in  such  a  way  as  to  cause  this  method  to  fail. 
Therefore  this  method  is  of  no  further  interest  in  this 
thesis. 

Many  of  these  methods  have  been  reviewed  by  Leese  and 
Boyd  (10).  They  recommend  an  integral  equation  technioue  as 
the  best  method  to  solve  the  M/M/1  queueing  problem  based  on 
computational  time  and  storage  considerations.  This  chapter 
draws  upon  their  review  and  computational  experience. 

The  six  more  useful  methods  are  now  discussed  in 
greater  detail.  In  the  first  four,  the  model  to  be  used  is 
ne  M/M/1  queue  with  equations  written  as 

P  (t)  =  —  o ( t ) P  (t)  +  p  (t)  ,  and  for  n>1 

0  0  1 

P  (t)  =  p(t)P  ft)  -  (o(t)+1)P  (t)  +  P  (t)  , 

n  n-1  n  n-H 


where 


o(t)  =  A  ( t ) / u  and  t  is  measured  in  units  of  ti . 
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Generating  Functions 

In  this  section,  transient  solutions  obtained  by 
application  of  generating  functions  are  discussed.  Using  the 
equations  for  Pn  (t)  as  defined  above,  the  following 
generating  function 

00 

F(x,t)  =  l  xnP  (t) 
n=0  n 

may  be  defined.  Partial  differentiation  with  respect  to  x 

and  t  and  elimination  of  P  (t)  terms  yields  the  Fokker- 

n 

Planck  Diffusion  equation 

X(1-X>  -  {1-X)2(1-PX)  ||  +  |£  +  Pd-X)2F  =  0. 

Solutions  to  this  equation  with  finite  difference 
techniques  are  frequently  difficult  to  obtain  due  to 
instabilities.  Other  techniques,  such  as  the  method  of 
characteristics  and  the  particle-in-cell  method  have  proven 
unsuccessful  in  various  degrees  when  p(t)  is  nonlinear  with 
respect  to  t. 

This  technique  has  only  thus  succeeded  in  transforminci 
a  difficult  problem  into  an  awkward  problem  for  the 
numerical  analyst.  Computationally,  partial  differential 
equations  also  require  large  amounts  of  core  storage,  making 
them  expensive  to  solve  as  well.  Also,  the  method  makes 
explicit  use  of  the  differential  equations  which  means  that 


lyfti  'tjtt u-'.-y^'y-yiy-V ’ KJJ  j.1  ti4|i I JHM&iB 
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for  each  model/  a  new  derivation  for  the  partial 
differential  equation  is  required. 

The  use  of  a  generalized  generating  function  of  the 

form 


F(x,t)  =  E  f ( x / t ) P  ( t ) 
n=0  n  n 

can  be  sometimes  used  to  advantage  when  a  particular  choice 
of  f  will  simplify  the  resulting  partial  differential 
equation.  For  example,  if  f  is  a  polynomial  in  x  with 
degree  n  and  P  (t)  is  as  above,  F(x,t)  will  sati.  :;y 


It  =P(X-1)F 

if  x  satisfies 


pfn+1  +  fn  +  £n- 1  =  (»x+1)fn  '  n  £  1 


with 


f  =  1  and  f  =  x. 
o  1 

In  this  case  the  second  order  partial  differential  eauation 
is  thus  simplified  to  a  first  order  system  with  an  auxiliary- 
algebraic  condition. 


-:  -1"  --■  .  i  m.i  i  ggcpeg TgfgP!&* g?  B!? 
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Another  means  of  solving  partial  differential  equations  is 
by  the  method  of  integral  equations.  This  method, 
sufficiently  important  in  its  own  right,  is  discussed  in  a 
separate  section  of  this  chapter  starting  on  page  13. 


Analytic  Continuation 

A  discussion  is  now  given  of  the  method  of  analytic 

continuation  of  the  function  P  (t)  in  the  interval  (t,t+h) 

n 

using  Taylor's  expansion. 

Knowing  P  (t)  for  some  t  allows  one  to  obtain 
n 

Pn(t)  for  all  n.  Repeated  differentiation  of  the 
equations  P^(t)  yields 


TTr|pn‘t>> 


Applying  Taylor's  formula,  the  expression  at  t+h  becomes 


P  (t+h) 
n 


h^ 

3T 


(j) 


dt 


T3T 


[p  (t)J  . 

n 


At  least  eight  terms  of  the  series  are  usually  required 
for  meaningful  accuracy.  We  need  only  note  that  if  p  is  a 
function  of  t,  the  chain  rule  of  differentiation  will  make 
calculation  of  the  derivatives  rather  tedious.  Although 
programs  for  doing  this  (FORMAC,  SNOBOL,  LISP,  etc.)  are 
available,  they  are  expensive  and  difficult  to  use. 
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The  use  of  this  technique  requires  an  error  analysis  of 
the  contribution  of  truncated  terms  of  the  series.  Usually, 
an  attempt  is  made  to  bound  the  error  with  an  expression 
involving  the  number  of  terms  before  truncation,  h,  and  the 
magnitude  of  higher  derivatives  of  P^  (t)  (1). 


Power  Series  Expansions 

This  technique  uses  the  Bessel  function  solution  which 
is  apparently  originally  due  to  Clark  (3,  15).  This  solution 
is  given  here  in  Chapter  three  on  page  24.  The  Bessel 
functions  and  exponential  terms  of  this  solution  are 
rearranged  such  that  Pn  (t)  is  expressed  as  a  power  series 
in  |0  of  the  form 


12 

This  method  works  well  if  p  can  be  approximated  by 
constant  values  over  short,  equally-spaced  intervals,  since 
the  C's  need  only  be  computed  once  and  stored  for  all 
N,m,n  and  t.  However,  a  large  amount  of  computing  time  and 
high-speed  storage  are  required.  Also,  the  method  is  not 
exact  and  requires  a  great  deal  of  effort  to  establish  the 
expressions  for  the  C's  if  the  form  of  the  equation 
changes.  The  problem  of  numerical  instability  in  evaluating 
the  terms  in  the  infinite  sums  must  be  constantly  considered 
by  the  analyst,  since  the  C's  can  exceed  unity. 

To  avoid  this  problem,  the  terms  in  the  power  series 
expansion  can  be  written  differently  leaving  the  Poisson 
part  unaltered.  Thus 


"  r(Pt)Vpt  m 
m=0  ml  N'n 


(t) 


where 


Fm 

N,n 


(t) 


S 


m+N-n 


(t)  + 


T 

m-n 


(t)  -T 

m-n- 

___ 

m 


1 


(t)  00 

l  S  (t). 
u=m+N+ 1  u 


In  this  form,  each  of  the  terms  f.]  and  F  represent 
probabilities  and  are  thus  less  than  unity.  Although  the 
possible  computational  difficulties  are  largely  avoided  with 
this  new  form,  obtaining  a  new  set  of  equations  for  each 
different  queueing  problem  encountered  would  be  tedious. 
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Integral  Equations 

In  this  section,  techniques  will  be  discussed  for 
solving  partial  differential  equations  such  as  the  form  of 
the  Fokker-Planck  equation  developed  in  the  section  under 
Generating  Functions.  A  solution  to  the  Fokker-Planck 
partial  differential  equation  may  be  formulated  in  terms  of 
a  Volterra-type  integral  equation?  however  it  is  not 
suitable  for  numerical  work. 

One  formulation,  depending  upon  a  change  of  variable, 
has  been  found  quite  suitable  for  solution  by  computer,  and 
in  fact,  is  probably  the  most  efficient  in  terms  of  speed 
and  storage  required.  This  is  the  method  of  Wragg  as 
implemented  by  Leese  and  Boyd  (10).  The  change  of  variable 
is 


q  (t)  =  l  P.  (t)  . 
n  j=n  ^ 


Then  by  substitution  into  the  equations  for  P  (t) 

n 


q  (t)  =  pq  (t)  -  (1+p)q  (t)  +  q  (t)  ,  n  >  1 
n  n-i  n  n-H  — 


with  q  (t)=  1  and  q  (0)  =  0  ,  n  >  1  . 

o  n  - 


By  defining 


u  =  f  P(y)dy  ,  v  =  [  ]^2  and  u  =  2  u(t-x)  ^ 2 

jx  r-x 
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and  writing  in  original  qn (t)  form,  we  obtain 


q„(t) 


=  |  tln-1  (w)p  (xjv1'"  S-In+1  (w)vn+1lexp(-u-t+x)G(x)dx. 


In  this  expression  I  (w)  is  a  Bessel  function  of  the 
First  Kind  and  G(x)  satisfies 


1  =  G (t)  +  f  K(t,x)G(x)dx  , 
*  o 


whore 


K(t,x)  =  [  -  v]  1^  (o))exp(-u-t+x)  . 

This  technique  still  involves  solution  of  Bessel 
functions,  complicated  by  the  fact  that  solution  for  G ( t ) 
requires  iteration.  The  above  integral  equation  is  far 
removed  from  the  original  model  in  terms  of  analytical 
transformations.  To  arrive  at  it,  the  form  and  assumptions 
of  the  original  queueing  problem  were  heavily  used.  Thus, 
altho  this  is  an  impressive  method  of  solving  for  the 
transient  conditions  in  the  M/M/1  model,  it  can  hardly  be 
relied  upon  as  a  technique  for  a  broad  class  of  problems. 

One  other  technique  using  partial  differential 
equations  for  solution  is  useful  in  large-scale  queueing 
systems.  This  method  is  due  to  Newell  (12)  who  applied  it  in 
highway  traffic  problems. 
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Define  F(x,t)  =  P[X(t)  <  x]  where  X(t)  is  the  queue 

length  at  time  t,  then  the  Fokk or- Planck  equation  may  be 

gp  (x  t) 

used  as  a  second-order  approximation  of  — -  where 


3F 

3t 


...  9F 
a(t)  3^ 


b(t)  3  2F 
2  3x2  * 


This  is  valid  for  largo  X,  a  larqe  number  of  arrivals  and 
departures  in  time  interval  t  with  only  a  small  change  of  X 
during  x  and  further,  that  the  number  of  arrivals  and 
departures  during  x  is  normally  distributed.  These 
conditions  are  often  not  unduly  restrictive  for  large  scale 
systems. 

Let  E  [A  (t)  3  be  the  number  of  arrivals  until  t  and 
E [D ( t) ]  be  the  number  of  departures  until  t. 

Then 


xa (t)  =  U [A ( t )  -  D ( t) ]  -  E (A ( t+x )  -  D(t+x)3 

=  E[D(t.+r)  -  D  ( t: )  +  A  (t+x)  -  A  (t)  ] 

and  similarly 

Xb (t)  =  VAR [D ( t+X )  -  D (t)  +  A (t+X)  -  A (t) 3 . 

By  defining  t*  =  ~  and  x*  =  —  then 

i  li 

0  0 
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and  F  (0)  =  (1-p.)F  (1)  ,  a«  1 

J+1  J  j  3 


which  can  be  applied  iteratively  to  solve 


3F .  (k)  3F.(k)  32F.(k) 

— I -  =  ( 1  -a . ) — - -  +  1 - 3 - 

3 j  3  3k  2  3k2 


Monte-Carlo  Simulation 

Monte-Carlo  simulation  is  a  powerful  tool  for  solving 
queueing  problems  in  the  sense  that  it  can  presumably  be 
applied  to  any  model.  In  cases  where  only  a  narrative 
description  of  the  problem  exists  or  where  differential 
equations  can  not  be  derived  describing  the  model,  Monte- 
Carlo  is  perhaps  the  only  quantitative  approach.  The 
relevant  question  then  is  how  efficient  is  Monte-Carlo  as 
compared  to  numerical  techniques?  Modern  computational 
equipment,  simulation  languages,  variance-reduction 
techniques  and  refined  cost  accounting  of  computer  programs 
have  made  Monte-Carlo  simulation  an  attractive  alternative 
for  an  analyst  who  may  not  wish  to  spend  a  great  deal  of  his 
time  attempting  to  derive  a  set  of  differential  equations  to 
describe  his  system.  Monte-Carlo  simulation  of  a  simple 
model  is  used  when  it  is  contemplated  that  the  model  will  be 
enriched  with  complexities  that  will  make  it  analytically 
intractable  at  a  later  time. 
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Monte-Carlo  simulation  is  a  reasonable  way  of 
estimating  the  steady-state  solution  of  a  system. 
Estimating  transient  behavior,  however,  is  much  mire 
difficult  and  costly  if  it  is  desired  to  estimate  P(t)  at 
many  discrete  values  of  t.  Since,  in  reality,  analysts  are 
generally  interested  in  obtaining  quantities  such  as  the 
waiting  time,  number  in  the  system,  and  the  number  of  units 
processed,  P(t)  is  not  required  and  the  desired  quantities 
may  oe  observed  directly  and  economically. 

A  simulation  of  a  M/M/1  queue  was  written  in  SIMSCRIPT- 
II  (.))  for  this  thesis  in  order  to  form  a  basis  for 
comparison  with  other  methods.  SIMSCRIPT-II  is  probably  one 
of  the  most  suitable  and  efficient  languages  available  for 
simulation  involving  transient  queueing  behavior.  The 
simulation  and  results  are  given  in  the  next  chapter. 


Runge-Kutta  Integration 

Direct  solution  of  the  differential  equations  by  Runge- 
Kutta  methods  is  not  only  simple  and  straightforward ,  but 
well  understood,  flexible  and  easy  to  apply  to  new  problems. 
The  numerical  analysis  and  behavior  of  Runge-Kutta  are  well- 
known.  The  method  is  the  most  stable  of  the  methods 
commonly  used  in  differential  equations  (Euler,  Milne,  Adams 
Bashforth,  etc.),  it  is  easy  to  use  and  produces  error  of 
order  h\  where  h  is  the  step  size.  This  author  has 
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written  a  subroutine  which  will  integrate  simultaneous 
differential  equations  of  arbitrary  size  and  order  and  has 
established  several  years  of  experience  with  it  on  wide 
classes  of  problems  encountered  in  differential  equations. 
The  Chapman-Kolmogorov  equations  are  a  particularly  "worry- 
free"  class  of  problems  solved  by  the  Runge-Kutta  method. 
Therefore,  it  seems  reasonable  to  attempt  the  development  of 
a  general  computer  program  to  solve  these  equations.  To  a 
large  extent,  this  lias  been  accomplished. 

The  author's  proqram  consists  of  unit  modules  with 

multiple  entry  points.  The  user  writes  a  subroutine  named 

INIT  with  only  the  minimal  amount  of  FORTRAN  coding 

necessary  to  describe  his  equations  and  problem 

peculiarities.  The  first  entrv  noint  is  used  for 

initialization.  The  second  entry  point  is  used  to  provide 

values  for  P  (t)  ,  given  P  (t)  and  t.  The  third  entrv 
n  n 

point  is  used  after  each  step  in  the  solution  to  allow  for 
computation  of  auxiliary  quantities  of  interest  such  as 
queue  statistics.  The  fourth  entry  point  is  used  at  the 
completion  of  the  solution  to  re- initialize  for  another 
problem  or  to  go  on  to  some  other  task.  These  entry  points 
will  be  discussed  in  greater  detail  in  Appendix  C.  The 
programming  of  subroutine  entry  points  is  discussed  in  (6) 
in  detail. 

The  Runge-Kutta  subroutine,  as  described  below,  has 
been  adapted  following  the  form  given  by  Nielsen  in  (14,  p. 
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178).  In  this  form#  which  is  one  of  many  alternative  forms 
in  the  literature,  it  is  possible  to  solve  a  set  of  n-th 
order  simultaneous  differential  equations  without  explicitly 
redefining  the  higher  order  derivatives  to  obtain  the  usual 
first  order  system. 

The  Runge-Kutta  subroutine  performs  four  steps  which 
involve  evaluation  of  the  differential  equations  to  move  the 
solution  for  Pfi(t)  from  t  to  t+At.  Initially  set 

AP.  =  P.  (t)  i=1 ,N. 

11 

Then  the  four  Runge-Kutta  steps  appear  as  below. 


R-K 

Step  No.  Operations  Performed 


k=1  t  =  t  + 

i  1 


•  At 

P.  *  P.+P.  -r  |  i*  ,N 

ll  1  1  2 


P  =  P  (t  )  ;  i=l,N  AP  «•  AP  +2P  ;  i=1,N 

ii  i  1  i  i  ii 


k=2  t  =  t  +  ^ 
2  1 


At 


P.+P.  ^  1  i=1,N 

l  ii  * 


P  =  P  (t  )  ;  i=1,N  AP  «•  AP  +2P  ;  i=1,N 

i2  i  2  i  i  i2 


k=3  t  =  t+  At 

3 


P.  =  At;  i=1,N 

is  i  12 
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P  =  P  (t  )  ;  i=1,N 
i  3  i  3 


AP  <■  -r-  (AP  +P  )  ;  i=1  ,N 
i  6  i  i s 


k=4  t  =  t  +  At.  P.  =  P.+AP.  ?  i=1,N 

4  14  1  1 


P  =P  (t  )  ?  i  =  1  ,N 
i1*  i  w 


t  =  t 


P  =  P  ; i~  1  ,  N  P  =  P  ; i=1  ,  N 

i  **  i  i  *♦ 


It  should  be  noted  that  only  three  workinq  storage 
vectors  of  length  N  need  be  provided  for  the  P^,P^  and 
P..  values.  These  are  provided  by  the  user's  subroutine 
and  may  be  used  for  other  storage  requirements  when  not  in 
use  by  tile  four-stop  Runqe-Kutta  computations. 

Rules  for  choosing  At  have  been  developed  to  allow  it 
to  change  as  the  solution  progresses,  depending  on  the 
magnitude  of  some  measure  of  truncation  error  (1). 
Generally,  these  rules  are  heuristic  or  highly  problem- 
dependent.  This  author  recommends  a  constant  At  initiallv 
and  then  re-running  the  oroblcm  at  a  few  selected  values  for 
At,  which  should  give  the  user  a  good  feel  for  the  trade-off 
between  computational  time  (cost)  versus  quadrature  error. 
Since  the  P^  values  must  sum  to  uni  tv  for  any  and  all  t, 
the  difference  between  this  sum  and  unity  makes  an  excellent 
measure  of  accumulated  error  at  each  step  of  the  solution. 

Further  work  on  this  problem  should  consider  the 
optimality  of  the  weights  of  the  P.-,. 


in  computing  the 
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AP^  above.  Currently, 

AP.  =At(aP.  +  bP.  +  cP.  +  dP.  )/w 
1  1  11  12  1 3 

where  a=1 ,b=2 ,c=2 ,d=1  and  w=a+b+c+d.  These  values  are 
optimal  if  the  function  being  integrated  is  a  parabola  on 
the  interval  from  t  to  t+At.  If  one  were  to  replace  this 
assumption  by  the  fact  that  the  P.  must  always  sum  to 
unity,  it  might  be  possible  to  derive  a  new  set  of  weights 
more  specifically  suited  to  this  kind  of  problem. 

The  next  chapter  will  present  a  detailed  comparison  of 
a  method  representative  of  those  discussed  in  this  chapter, 
the  Runge-Kutta  method  and  Monte-Carlo  simulation. 
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CHAPTER  3 

A  COMPARISON  OF  THREE  METHODS  FOR  THE  M/M/1  QUEUE 

In  this  chapter  a  comparison  is  made  between  a 
representative  analytic  solution,  the  Runge-Kutta  solution 
of  a  single-channel,  single-server  queue,  and  the  Monte- 
Carlo  method.  The  analytic  solution  for  this  model  is 
presented  in  Saaty  (15)  and  is  considered  representative  in 


complexity 

and  computational 

difficulty  of 

all 

the 

analytical 

methods  discussed 

in  the  previous  chapter. 

The 

comparison 

of  methods  will  consist  of  computing 

P(t) 

from 

t=0  until 

• 

steady  state  is 

achieved,  i.e., 

t=t* 

where 

Analytic  Bessel  Function  Method 
The  time-dependent  solution  of  the  single-channel, 
single-server  queue  is  derived  (apparently  by  Clarke  (3))  in 
terms  of  modified  Bessel  functions  of  the  First  Kind. 


24 


where  i  is  the  initial  number  of  items  in  the  system.  A 
computer  program  was  written  in  FORTRAN  using  the  Bessel 
function  subroutine  BESI  from  the  IBM  Scientific  Subroutine 
Package  (SSP)  (7).  The  listing  of  this  program  with  sample 
output  is  included  in  Appendix  B.  The  size  of  the  load 
module  was  29,856  bytes.  Time  t  was  varied  from  0  to  20  in 
steps  of  0.1  and  n  varied  from  0  to  20  with  i=0.  The 
parameters  A  and  u  were  chosen  as  0.45  and  0.5,  res¬ 
pectively.  The  equation  for  P  (t)  above  was  solved  for  201 

n 

values  of  t  and  21  values  of  n  for  a  total  of  4221 

times.  The  amount  of  central  processor  unit  (CPU)  time  for 

the  IBM  360/65  was  98.64  seconds.  Approximately  four  hours 

were  required  to  write  the  program  and  several  trials  were 

needed  to  debug  it.  The  infinite  series  term  of  the 

equation  was  truncated  when  either  k=n+50  or  the  k-th  term 

was  less  than  10 This  rule  produced  values  of  P  (t) 

n 

which  were  no  more  than  about  IX  in  error. 
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Runge-Kutta  Method 

The  single-channel,  single-server  equations  were 
integrated  using  the  Runge-Kutta  method.  For  completeness 
they  are  listed  here  again  as 

P  (t)  =  -AP  (t)  +  UP  (t) 
o  o  i 

P  (t)  =  AP  ,  (t)  -  (A+u ) P  (t)  +  UP  ,  (t)  for  1  <  n  <  N 
n  n-i  n  n+1  ~ 

P  (t)  =  AP^ft)  -  (A+u) P^ (t)  with  PQ  (0)  =  1. 

The  same  values  for  t,A  and  u  were  used  with  N=20.  The 
time  step  At  was  chosen  as  0.1.  Since  the  Runge-Kutta 
method  requires  four  evaluations  of  the  equations  per  step, 
the  above  equa cions  were  evaluated  a  total  of  201  x  4  =  804 
times  for  each  value  of  n.  This  is  804  x  21  =  16,884 
equations  evaluated  for  the  run.  The  Runge-Kutta  method 
requires  storage  of  the  entire  P(t)  and  P(t)  vectors  plus 
three  more  working  vectors  of  the  same  length.  The  size  of 
the  load  module  for  this  program  was  25,960  bytes.  The 
amount  of  CPU  time  for  the  solution  was  7.42  seconds.  The 
solution  error  grew  to  a  maximum  of  0.005X  with  this  step 
size.  Only  the  differential  equations  and  initial  conditions 
had  to  be  programmed  as  a  subroutine  named  INIT  since  the 
main  program  is  completely  general.  This  consumed  about  1/2 
hour  and  did  not  need  to  be  debugged. 
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The  following  table  illustrates  how  CPU  time  varies  as 
a  function  of  the  maximum  number  in  the  system  N  for  the 
Bessel  function  solution  and  for  Runge-Kutta  solution.  Times 
are  in  CPU  seconds. 

Table  1 

Comparative  CPU  Times 


N 

Bessel  Functions 

Runge-Kutta 

10 

59.0 

5.0 

20 

100.9 

7.4 

30 

149.9 

9.6 

40 

205.2 

11.3 

50 

290.6 

12.9 

75 

677.9* 

16.4 

100 

1428.5* 

21.8 

The  asterisk  (*)  notation  above  indicates  these  figures  are 
linear  extrapolations  based  upon  runs  over  a  shorter  range 
of  t,  since  there  is  a  uniform  correspondence  between  CPU 
time  for  the  run  and  t. 

When  plotted,  the  times  for  the  Bessel  function 
solution  vary  as  power  function  with  N  (Pigure  1).  The 
Runge-Kutta- solutions  vary  linearly  with  N  (Figure  2). 


FIGURE  1 


CPU  Tim 
Vtrtui  Nunber  in 
the  Syetea  M  for 
Seieel  Function  Method 
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Estimates  of  required  CPU  time  T  may  be  obtained  from  the 
empirical  relations 

T  =  46  exp(0.0355N) 
for  Bessel  functions,  and 

T  =  4.2  +  0.172N 

for  Runge-Kutta.  These  were  derived  from  the  plots  of  the 
above  CPU  times.  The  cost  of  the  additional  core  storage 
required  to  use  the  Runge-Kutta  method  beyond  the  break-even 
point  of  N=500  is  insignificant  in  comparison  to  the 
increased  CPU  time  cost  required  by  the  solution  by  Bessel 
functions. 

The  FORTRAN  subroutine  INIT  for  the  M/M/1  queue 
discussed  in  this  chapter  is  presented  in  Appendix  C  along 
with  sample  output  used  in  the  comparison. 

The  following  table  is  a  summary  of  benefits  realized 
by  using  the  Runge-Kutta  method  over  the  analytic  solution 
involving  Bessel  functions.  Each  benefit  more  or  less 
applies  to  every  analytical  technique  discussed  in  the 
previous  chapter. 
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Table  2 

Summary  of  Benefits 

Runge-Kutta  Bessel  Functions 


analytical 

No  further  analysis 
required  beyond  derivation 
of  equations. 


difficulty 

Considerable  analysis  to 
get  analytic  solution  to 
equations.  Analytic 
expression  may  not  exist  - 
may  need  approximations. 


programming 

Easy  to  program  for  Sophisticated  programming 

computer  solution  -  little  expertise  required  to 

debugging  required.  prevent  overflows, 

loss-of-significance,  etc. 
Difficult  to  debug. 


numerical 

Solution  can  be  trusted, 
method  behaves  in  a  pre¬ 
dictable  way,  precision 
can  be  arbitrarily  chosen 
by  user. 


stability 

Numerical  instabilities 
leave  trustworthiness  of 
solution  in  doubt,  very 
little  flexibility  in 
choosing  desired  precision. 
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cost  of  accuracy 


Cost  of  achieving  a  given 
level  of  precision  is 
related  linearly  to  step 
size  over  wide  range  of 
values. 


Cost  of  achieving  a  given 
level  of  precision  usually 
is  related  ncn-linearly  in 
a  power  function  over  a 
narrow  range  of  precision. 


analysts'  time 


Time  to  obtain  solution, 
given  a  model,  is  short. 


Time  to  obtain  solution, 
given  a  model,  may  be  very 
long. 


clarity  of  methods 


Method  and  solution 
easily  understood  by 


Method  chosen  and  imple¬ 
mentation  difficult  for 


others . 


others  to  understand. 


cost  of  solutions 


Cost  of  solution  increases 


Cost  of  solution  increases 


linearly  with  N  ,  the 
number  of  equations 
solved. 


as  a  power  function  of  N. 


size  of  programs 


Load  module  of  equal 
size  occurs  at  approxi- 


Load  module  size  does  not 
change  with  N;  greater  for 


mately  N=500. 
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N  <  500  than  the  Runge- 
Kutta  load  module. 


Monte-Carlo  Simulation 

A  simulation  of  the  M/M/1  model  was  written  in 
SIMSCRIPT-II  (9).  The  listing  and  sample  output  are 
contained  in  Appendix  A.  For  a  sample  size  of  1000 
replications,  consisting  of  10  runs  of  100  replications 
each,  the  95%  confidence  interval  on  the  expected  number  in 
the  system  n  was  obtained  as 
1.66  <  n  <  2.32. 

Although  this  is  quite  a  large  variability,  it:  is  acceptable 
to  achieve  this  level  in  complex  simulations.  To  make  a 
comparison  with  the  other  two  methods  on  the  basis  of 
computer  resources,  the  load  module  size  was  48,000  bytes 
and  the  total  CPU  time  was  28  seconds.  This  compares 
favorably  with  Runge-Kutta  and  shows  the  efficiency  of  the 
SIMSCRIPT-II  system.  The  CPU  time,  of  course,  would  depend 
linearly  on  the  number  of  replications.  The  time  to  program 
and  debug  the  simulation  was  slightly  less  than  that 
required  for  the  Bessel  function  program. 

Although  not  strictly  a  numerical  technique,  the 
simulation  was  included  to  demonstrate  its  relationship  with 
the  previous  two  methods  because  simulation  is  usually 


regarded  as  an  expensive  and  inefficient  method  to 
such  problems. 


solve 
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CHAPTER  4 

EXAMPLES  OF  SOLUTIONS  TO  OTHER  QUEUEING  MODELS 

In  this  chapter  two  other  applications  of  the  Runge- 
Kutta  method  will  be  presented.  The  purpose  is  to  show  how 
simple  and  straightforward  the  preparation  for  solution  of 
different  problems  can  be  when  using  the  Runge-Kutta  control 
program  in  Appendix  C.  The  first  is  the  solution  of  the 
premptive-resume  priority  queue  and  the  second  is  the 
solution  of  the  state-dependent  queue.  These  were  chosen  as 
further  examples  of  solution  by  the  Runge-Kutta  integration 
because  they  are  representative  of  commonly-used  models. 


Premptive- Resume  Queue 

The  premptive-resume  queueing  model  is  a  slight 
variation  of  the  M/M/1,  In  it,  the  possibility  of  a 
different  kind  of  priority  arrival  is  allowed  in  the  system 
with  a  different  arrival  rate  than  the  normal  arrivals. 
When  a  priority  arrival  occurs,  service  is  suspended  if 
service  is  being  performed  on  a  normal  service  item,  and  it 
is  set  aside  so  that  service  may  begin  on  the  priority 
arrival.  Service  proceeds  at  a  different  rate  for  the 
priority  item  than  for  normal  items.  The  server  continues 
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to  serve  from  the  priority  queue  until  empty,  at  which  time 
the  normal  item  that  had  been  set  aside  will  return  to 
service  and  continue  where  it  left  off. 

This  is  an  excellent  model  for  describing  how  a  M/M/1 
model  would  behave  if  the  server  were  subject  to  breakdown 
and  repairs  were  made  to  restore  it  to  service. 

In  the  model  that  follows  X  and  u  are  the  arrival 
and  departure  rates,  respectively,  for  the  normal  system  and 
Xp  and  hp  are  parameters  for  the  priority  system.  The 
equations  are  written  for  a  finite  number  M  in  the  priority 
system  and  N  for  the  normal  system.  Moreover  the  subscript 
i  stands  for  the  number  in  the  priority  system  and  j 
stands  for  the  number  in  the  normal  system. 
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pi.i(t)  ■  -(x+V“p)pi,j<t>+Vi-.,j<t)  +  + 

P0j0(t)  .  -a+Ap)P0>0(t)  +  gpPlf0(t)  +  uP|)>i(t) 

PCjJ(t,  .  -(A+Ap+u)Po  j(t,  +  AP0fJ_l(t)  +  gP0ij+1(t)  +  UpPij(t) 

pi,0(t)  ■  -(>+V"p)pi,o(t)  +  +  g„P1+1>0(t) 

■  -«+V+“p)PM,N(t>  +  Vh-u«M  +  XPM.l.-,(t> 

po,N(t)  ■  -<»+y»>pM<t)  +  »P0iH.,(t)  +  MpP11((t> 

PM,0(t)  "  '(>+Xp%)PM,0(t)  +  VM-l,0(t> 

'm.J  (t>  =  -a+V"p>PM,J<t>  +  XpPM-],j(C)  + 


Subroutine  INIT  for  use  with  the  program  in  Appendix  C  and 

sample  output  follow,  along  with  a  plot  of  some  auxiliary  queue 
statistics  in  Figure  3. 
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State-Dependent  Queue 

The  model  discussed  here,  is  a  generalized,  sinqle- 
channel,  state-dependent  queue  with  s  servers  developed  by 
Conway  and  Hillier  (4).  The  differential  equations,  in 
terms  of  the  independent  variable  x  =  pt  ,  are 

P  (t)  =  n  P  (t)  -p  (r)  P  (t) 

0  11  oo 

p  (t)  =  p  .(iJP  .  (t  )  -  [p  (x)+n  ]p  (x)+n  p  .(x)  , 
n  n-1  n-1  n  n  n  n+1  n+1 

for  1  <  n  <  N 

p  (x)  =  P  ( T ) P  (x)-[p  (x)+n  )  P  (x)  , 

N  n-1  N-1  N  N  N 

where  N  is  the  maximum  number  of  items  in  the  system, 


p  (x)  =  X  (x)/p  , 

n  n 


n 


n 


=  V  /v 

n 


n  <  s 


i 


n  >  s 


A  (t)  ,  n  <  s 


A  (t)  = 
n 


(HiT)bA{t)  ,  n  >  s. 
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This  model  becomes  the  standard  M/M/1  queue  when  the 
parameters  b  and  c  are  zero  and  s=1.  This  is  a  useful 
model  for  working  with  an  existing  system  for  which 
empirical  data  are  available.  The  time  dependence  of 
arrivals  and  departures  can  be  handled  as  well  as  a  wide 
variety  of  system  dependence  through  the  choice  of  b,  c  and 
s . 

Here  p  =  A/y  is  called  the  loading  factor.  A  steady- 
state  solution  can  not  be  obtained  when  P>1.  However,  it  is 
often  required  to  study  systems  for  which  p  _>  1  during  part 
of  the  time  of  interest.  Figure  4  shows  how  a  single  channel 
dual  server  queue  model  behaves  with  p  =  1.5  from  t  =  n  until 
t  =  8  with  the  system  initially  empty. 

Subroutine  INIT  for  this  model  and  sample  output  are 
included  next.  Note  that  T  in  the  output  refers  to  t 
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CHAPTER  5 
SUMMARY 

Methods  for  transient  solutions  of  queueinq  problems 
have  been  explored  and  evaluated  for  their  usefulness  in 
numerical  work.  The  so-called  closed-form  analytic 
solutions  are  the  least-preferred  forms  for  numerical  work 
based  on  difficulty  of  analysis,  computer  programming  and 
cost.  Monte-Carlo  simulation  is  easier  and  cheaper  to  apply 
than  these  analytic  solutions  but  should  not  be  used  when 
the  queueing  model  can  be  expressed  in  the  form  of 
differential-difference  equations,  In  this  case,  the  use  of 
Runge-Kutta  integration  is  preferred. 

The  main  control  program  in  Appendix  C  has  been  shown 
to  be  easy  to  use,  capable  of  high  precision,  low  in  cost 
and  useful  in  a  wide  range  of  queueing  models.  It  is 
offered  as  a  much-needed  tool  to  researchers  and 
practitioners  in  queueing  theory  with  the  hope  that  the  use 
of  transient  analyses  in  queueing  theory  will  be  more  easily 
attainable.  To  this  extent,  existing  results  will  be  made 
more  meaningful  (2). 

Two  areas  of  further  research  have  become  apparent 
during  the  preparation  of  this  thesis.  One  concerns  the 
incorporation  of  the  Runge-Kutta  queueing  problem  solver 
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described  in  Appendix  C  in  continuous  ootimization 
applications.  The  other  concerns  numerical  analysis  on  the 
Runge-Kutta  algorithm  as  employed  here. 

The  ability  to  obtain  fast  and  inexpensive  transient 
solutions  to  queueing  problems  suggests  the  application  of 
variational  methods  to  compute  optimal  functional  solutions 
rather  than  simple  scalar  solutions.  An  example  might  be  a 
problem  of  finding  the  optimal  .service  rate  of  a  system 
subject  to  a  given  fluctuating  rate  of  arrivals  and  to 
constraints  on  the  maximum  number  in  the  system.  The 
optimality  conditions  might  include  a  cost  function  relating 
how  cost  varies  as  a  function  of  service  rate  or  capacity, 
number  in  the  system  and  the  loss  of  arrivals  when  the 
system  is  full.  The  optimizing  equation,  representing  the 
integral  of  the  variational  of  the  optimality  conditions 
would  then  have  the  solution  of  the  queueing  model  imbedded 
within  its  solution.  A  new  main  program  to  solve  the 
queueing  problem  as  a  part  of  the  optimization  solution 
would  be  required  in  addition  to  what  is  presented  in  this 
thesis.  Numerous  optimization  algorithms  are  available  for 
ready  incorporation  in  such  a  program. 

Additional  work  in  the  numerical  analysis  of  the  Runge- 
Kutta  algorithm  was  suggested  in  Chapter  two.  The  Runge- 
Kutta  algorithm  can  be  derived  in  many  forms  depending  upon 
the  desired  error  in  a  single  step  and  assumptions  regarding 
the  behavior  of  the  equations  being  integrated.  The  form 
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employed  in  this  thesis  is  one  commonly-used  on  a  qreat 
variety  of  problems.  It  can  be  made  to  work  in  most 
difficult  cases  by  decreasing  the  step  size  until  stability 
is  eventually  achieved.  Generally,  when  this  occurs,  the 
analyst  realizes  the  algorithm  may  not  be  aonrooriate  and 
seeks  another  more  suited  to  the  problem. 

On  the  other  hand,  it  has  been  discovered  that  the 
Runge-Kutta  equations  will  solve  the  queueing  models  in  this 
thesis  accurately  with  very  large  step  sizes.  This  would 
seem  to  indicate  that  the  method  is  more  accurate  than  it 
need  oe.  The  differential  equations  in  queue ing  are  a 
specialized  type  of  problem.  They  are  usually  first-order, 
nighly-coupled  and  yet  stable  under  drastically-changing 
conditions  of  the  rate  parameters.  Future  work  should 
attempt  to  rederive  the  Runge-Kutta  formula  for  the  oueuoina 
equations.  It  may  be  possible  to  obtain  a  more  economical 
form  utilizing  tiie  characteristics  of  aueuemc  models  and 
with  fewer  evaluations  of  the  differential  equations  at  each 


step. 
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APPENDIX  A 

LISTING  OF  MONTE-CARLO  SIMULATION  PROGRAM 

Tiie  computer  program  listed  here  is  a  simulation  of  a 
M/M/1  queue  with  the  same  rate  parameters  and  maximum  number 
of  items  in  the  system  N  as  used  in  the  comparison  in 
Ciiapter  three.  The  operation  of  the  simulation  is 
straightforward  and  largely  self-explanatory  due  to  the 
uniqueness  of  the  SIMSCRIPT-II  language.  The  elements  of 
tiie  simulation  include  a  QUEUE  which  is  FIFO  by  default, 
ARRIVAL  events  which  cause  ITEMS  to  be  created  and  SERVICE 
events.  When  an  ARRIVAL  event  occurs,  an  ITEM  is  created, 
its  time  of  arrival  is  saved  and  it  is  filed  in  the  QUEUE. 
The  time  between  arrivals  is  an  exponentially-distributed 
random  variable  computed  by  the  EXPONENTIAL.?  function. 
ARRIVALS  reschedule  themselves.  A  SERVICE  event  is 
triggered  whenever  an  ITEM  joins  th'  pty  QUEUE  and/or  when 
an  ITEM  completes  SERVICE  a^d  the  QUEUE  is  not  empty.  The 
SERVICE  completion  event  time  is  also  scheduled  using 
exponentially-distribu'  1  random  variables. 

At  pre-specified  times  during  the  t^me  the  simulation 
is  allowed  to  run,  events  named  LOOK  occur.  This  event 
records  the  instantaneous  number  in  the  system  in  a 
histogram  named  P.  The  simulation  also  computes  averaae 
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waiting  time  in  the  system  and  percent  utilization  of  the 
server. 

The  simulation  is  replicated  1000  times  while  the 
histogram  is  accumulated.  The  histogram  cells  are  then 
divided  by  1000  to  give  estimates  of  the  probability  of  the 
number  in  the  system  at  a  particular  time  and  printed.  The 
sample  output  shows  this  for  times  of  5,  10/  15  and  20. 
These  times  correspond  to  times  printed  in  the  other  two 
methods  discussed  in  Chapter  three.  The  results  are  in  fair 
agreement  with  the  Runge-Kutta  solution  in  Chapter  three. 
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APPENDIX  B 

LISTING  OF  BESSEL  FUNCTION  PROG  PAM 

The  following  is  a  FORTRAN  listing  and  sample  output 
of  a  program  written  specifically  to  obtain  the  numerical 
solution  of  the  Bessel  function  transient  analytical 
expression  for  the  M/M/1  queue  (15). 
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APPENDIX  C 
RUNGE-KUTTA  PROGRAM 

The  Chapman-Kolmogorov  equations  representing  a 
queueing  model  are  readily  solvable  using  the  Runge-Kutta 
method  of  solving  differential  equations.  Experience  has 
shown,  however,  that  a  technique  of  standardization  of  these 
equations  and  the  form  of  the  Runge-Kutta  algorithm  to  solve 
them  must  be  sought  in  order  to  obtain  generality, 
convenience  and  efficiency.  This  section  describes  a 
programming  package  which  will  guide  its  user  in  obtaining  a 
systematic  time-dependent  solution  of  the  differential 
equations  of  queueing  models. 

It  is  assumed  that  the  user  of  this  programming  package 
is  familiar  with  the  Chapman- Kolmogorov  form  of  queueing 
equations  and  has  a  particular  problem  to  be  solved. 
Furthermore,  it  is  assumed  that  the  user  is  familiar  with 
the  FORTRAN-IV  programming  language. 

The  user  should  be  acquainted  enough  with  the  Runge- 
Kutta  algorithm  to  realize  that  it  is  a  self-starting  method 
which  requires  four  evaluations  of  the  differential 
equations  per  integration  step  h  in  the  independent  variable 
t.  Preliminary  evaluations  are  made  twice  with  the 
independent ' variable  at  fc+h/2,  one  at  t+h  and  then  one  final 
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evaluation  with  correct,  updated  solutions  at  t+h.  In  the 
program,  these  steps  are  referred  to  by  1,  2,  3  and  4.  The 
variable  k  will  also  be  used  to  stand  for  these  numbers. 
The  term  self-starting  above  means  that  only  values  of  the 
solution  variables  and  differentials  at  time  t  are  needed  to 
obtain  a  new  solution  at  time  t+h. 

With  these  preliminaries  in  mind,  it  is  n ow  possible  to 
provide  a  concise  summary  of  the  basic  scheme  which  will  be 
presented.  To  obtain  numerical  solutions  of  differential 
equations  using  the  Runge-Kutta  method  requires  one  to 
specify  the  equations  in  FORTRAN  notation,  specify  the 
initial  conditions,  specify  values  for  constants  and 
parameters  as  required  by  the  particular  problem  and  to 
specify  the  range  of  the  independent  variable  t  and  the  step 
size  h.  These  conditions  must  be  met  by  the  user  for  any 
problem. 

Normally,  in  addition  to  specifying  conditions  relating 
directly  to  his  problem,  the  user  must  implement  the 
integration  scheme,  define  some  type  of  data  structure  for 
intermediate  or  working  storage  and  construct  an  output 
format  and  control  over  print  frequency. 

The  programming  package  presented  here  standardizes  the 
method  of  defining  the  necessary  conditions  for  any  problem 
mentioned  formerly  and  thus  eliminates  the  necessity  to 
plan,  program  and  debug  the  latter  steps  relating  to 
implementation  detail.  This  is  accomplished  by  requiring 
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the  user  to  construct  a  subroutine  named  INIT  with  FORTRAN 
Entry  points  named  PRIME,  NEWVAL  and  WRAPUP.  This 
subroutine  will  describe  the  user  problem  in  standard  form 
to  the  programming  package. 

Briefly,  initial  entry  point  INIT  will  be  used  by  the 
user  to  specify  initial  conditions  of  the  problem,  to 
specify  the  range  of  the  independent  vari:j/»le  t  and  the  step 
size.  These  data  will  be  assigned  to  vav" ables  listed  in  a 
specified  order  in  a  C0MM0N  statement..  Entry  point  PRIME 
will  serve  to  provide  computed  values  o*  the  differentials, 
given  values  of  the  solution  variable*  during  operation  of 
the  generalized  Runge-Kutta  subroutine  Through  entry  point 
NEWVAL,  the  package  makes  the  solutio'  variables  and  updated 
differentials  available  after  each  s^f.  h  in  t.  Entry  point 
WRAPUP  is  used  after  the  problem  is  ^clved  to  allow  the  user 
to  decide  whether  to  start  a  new  problem  or  terminate  the 
computer  run. 

Table  three  Program  Logic  presented  below  shows  how  the 
user  subroutine  INIT  interfaces  with  the  programming 
package,  which  consists  of  a  Main  program  and  the  Runge- 
Kutta  subroutine.  It  is  not  necessary  for  a  user  to 
understand  the  programming  package  itself,  but  only  to 
understand  the  responsibilities  required  of  the  INIT 
subroutine  at  each  of  the  points  of  interface. 

Following  the  Program  Logic,  the  next  part  of  this 
section  will  consist  of  a  detailed  example  of  how  the 
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equations  for  the  M/M/1  queue  were  solved  using  this 
systematic  programming  package. 
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Solution  of  the  M/M/1  Queue 
The  equations  for  the  M/M/1  queue  are 
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£  (t)  =  -AP  (t)  +  pP  (t)  , 
o  o  1 

Pn{t)  =  AP^Ct)  "  (X+u)P  (t)  +  yp  +1(t)»  for  1  <  n  <  N, 

PM(t)  =  APm  .(t)  -  (A+p)  P  (t) 

N  N-1  N 

where  N  is  the  maximum  number  of  items  in  the  system. 

In  this  example  P  (0)=1  and  P  (0)=0  for  n=1?N  are  the 

o  n 

initial  conditions.  The  parameters  are  A=0.45,  p=0.5  and 

N=20 .  The  variable  t  will  run  from  0  to  20  with  Runge-Kutta 
step  size  of  0.1. 

Although  theoretically  there  is  no  reason  to  limit  the 
maximum  number  in  the  system,  in  actual  computation  an  upper 
bound  must  be  chosen.  One  can  either  truncate  negligible 
terms  or  modify  the  equations  to  take  care  of  the  finite 
waiting  room,  as  was  done  above.  The  latter  choice  is 
recommended  because  it  is  then  possible  to  utilize  the 
auxiliary  condition  * 

i 

00 

2  P  (t)  =  1 
n=C  n 

to  compute  a  measure  of  solution  error  e(t)  from 
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N  A 

e (t)  ~  1  -  Z  P  (t)  . 
n=0  n 

A 

The  P  (t)  notation  is  used  to  distinguish  between  the 
n 

notation  used  in  the  model  and  the  actual  computed  values. 

The  usual  practice  is  to  assign  initial  conditions  at 
t=0  and  obtain  the  solution  of  the  equations  at  selected 
positive  values  of  t.  The  programming  package  will  print 
the  solution  after  every  multiple  of  NPRNT  Runge-Kutta  steps 
chosen  by  the  user. 


Subroutine  INIT 

At  this  point,  we  have  a  statement  of  the  necessary 
conditions  of  the  sample  problem  and  are  ready  to  write  a 
subroutine  INIT  for  solving  it.  The  following  will  be  a 
step-by-step  procedure  for  writing  the  FORTRAN  subroutine 
INIT.  Each  FORTRAN  statement  will  be  followed  by  an 
explanation  of  it  in  relation  to  the  problem  to  be  solved. 
Statements  marked  by  /  indicate  the  statement  is  always 
required  by  the  programming  package  interface  for  any 
problem,  except  for  the  characters  underscored,  which  are  to 
be  provided  by  the  user  as  required  by  the  problem. 
Statements  not  so  indicated  are  user-specified  to  pertain  to 
the  problem  at  hand. 

/  1)  SUBR0UTINE  INIT 

/  2 )  C0MM0N  NAME ( 2  0 ) , T0 , TM AX , DT , NP  RUT , NEQ , 
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P  (2 1 ,2)  ,W<63) 

This  statement  is  used  as  a  common  reference  to  locate 
data  for  the  problem.  Each  array  except  W  and  each  variable 
in  C0MM0N  must  be  assigned  values  as  described  below. 

3)  REAL  LAMBDA, MU 

Indicates  X  and  y  ,  respectively,  in  the  equations  will 
have  these  names  and  will  be  floating  point  reals. 

✓  4)  READ  (5,1,  END=  1 0 )  NAME ,  LAMBDA ,  MU ,  TMAX ,  DT 

/  5)  1  F0RMAT  (2QA4/4F1 0.0) 

Reads  a  problem  identification  card  from  the  card 
reader  and  stores  in  array  NAME .  Reads  the  parameters  X  and 
y  ,  the  upper  limit  of  time  t  for  the  equations  TMAX  and  the 
step  size  for  Runge-Kutta  integration  h,  named  DT  in  the 
program. 

6)  PRINT  2, -NAME, LAMBDA, MU, TMAX, DT 

7)  2  F0RMAT  ( 1  Hi  20A4/1H0  9X6HLAMBDA  13X2HMU 

1 1X4HTMAX  1 3X2IIDT/1H0  4F15.4) 

Prints  out  the  input  data. 

✓  8)  T0=O . 0 

Assign  the  initial  value  of  time  t,  named  T  in  the 
program.  This  could  have  also  received  its  value  in  the 
READ  statement  as  well.  The  problem  is  to  be  integrated 
from  T0  to  TMAX  in  steps  of  DT. 

/  9)  NPRNT=5 

This  is  the  number  of  DT  steps  in  time  t  between 
printouts  of  the  solution  by  the  package  Main  program. 


/  10) 
The 


NEQ=2 1 
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problem  has  21  equations  to  be  solved.  Since 
FORTRAN  does  not  allow  a  subscript  of  zero,  all  subscripts 
of  the  problem  must  be  increased  by  one  for  the  program. 


11) 

NM0=NEQ- 1 

12) 

NP0=NEQ+1 

13) 

C0NS=LAMBDA+MU 

Constants  used  by  this  problem  to  control  loops  and 
prevent  needless  computation  in  later  entry  point  sections 
of  the  subroutine. 

/  14)  D0  3  1=2, NEQ 

/  15)  3  P(I,1)=0.0 

/  16)  P(1,1)=1.0 

Initialization  of  the  P-array.  A  later  section  titled 
Dimensioning  the  Problem  explains  the  reason  for  an 
additional  subscript  for  the  P-array. 

17)  LINE=55 

18)  IPRNT=0 

Auxiliary  counters  used  by  the  problem  program  for 
additional  output  control. 

/  19)  RETURN 

Returns  to  the  Main  program  to  complete  problem 
initialization. 
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PRIME  Entry  Point 

This  is  the  entry  point  called  from  the  Runge-Kutta 
suoroutine  during  the  four  step  estimation  of  the  solution 
vector.  The  estimates  are  in  an  array  named  PT  which  is 
given  to  this  entry  point  along  with  time  T  and  the  Runge- 
Kutta  step  number  K.  The  user  must  use  the  values  of 
PT(*,1)  in  the  differential  equations  which  are  stored  in 
PT(*,2).  The  *  notation  here  stands  for  the  subscripting 
pertaining  to  the  users*  problem,  explained  below  in  the 
section  titled  Dimensioning  the  Problem. 

J  20)  ENTRY  PRIME (T,P? , K) 

/  21)  DIMENSI0N  PT(21,2) 

The  PT-array  must  be  dimensioned  exactly  as  the  P-array 
was  in  the  C0MM0N  statement  above. 

/  22)  PT  (1_,  2)  =RHS  (  iiU*PT  (2,1)  -LAM3DA*PT  (1,1)  ) 

First  user  problem  equation  solved  and  stored  in  upper 
half  or  P-array.  Enclosing  the  expression  in  parentheses  as 
an  argument  to  the  RHS  fur  ction  prevents  machine  underflow 
interrupts. 

/  23)  D0  4  1=2 , NM0 

Loop  counter  for  the  general  equation  of  user  problem. 

/  24)  4  PT(I,2)=RHS{  LAMBDA*PT (1-1,1) -C0NS *PT (I, 1 ) 

+MU*PT  (1+1,1)  ) 

The  general  equation  of  the  user  problem  is  computed 
and  stored  in  the  upper  half  of  the  PT-array. 

/  25)  PT(NEQ,2)=RHS (  LAMBDA*PT (NM0, 1 ) -C0NS*PT (NEQ, 1 ) 
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The  final  equation  for  the  user  problem  is  computed  and 
stored  in  the  upper  half  of  the  PT-array. 

/  26)  RETURN 

Return  control  to  the  Runge-Kutta  subroutine. 

NEWVAL  Entry  Point 

This  entry  point  is  called  after  each  step  in  the 
solution,  i.e. ,  after  each  time  that  a  new  solution  vector 
is  computed  and  T+T+DT.  The  user  is  given  the  new  value  of 
time  T  and  the  solution  vector  and  updated  derivatives  in 
the  PT-array. 

/  27)  ENTRY  NEWVAL (T,PT) 

In  this  problem  it  is  desired  to  compute  and  print  the 
solution  error  after  every  NPRNT  steps  of  DT  in  T  and  label 
the  output  at  the  top  of  each  page. 

28)  CALL  C0NTRL ( IPRN1 , NPRNT ,  &  1 1 ) 

C0NTRL  is  part  of  the  programming  package.  Using  the 
free  counter  variable  IPRNT  and  counter  limit  NPRNT,  it 
returns  to  statement  11  whenever  IPRNT  is  less  then  NPRNT 
after  incrementing  IPRNT  by  one.  It  returns  to  the  next 
statement  if  IPRNT  equals  NPRNT  after  setting  IPRNT  to  zero. 

29)  CALL  C0NTRL (LINE, 55 , 65) 

This  counter  switch,  used  for  page  control,  is  similar 
to  the  above.  It  executes  the  next  statement  after  every  55 
passes . 

30)  PRINT  7, NAME 
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Print  column  headings  at  top  of  each  page,  that  is, 
after  every  55  passes  through  this  section. 

31)  7  F0RMAT  (1H1  20A4/1H0  1 4X  1HT  6X 

14HERR0R,  PERCENT) 

32)  5  SUM=0 .0 

Initialize  for  computing  the  solution  error. 

33)  D0  8  1=1 ,NEQ 

Loop  fo^  computing  sum  of  solution  vector. 

34)  8  SUM=SUf!+PT  (NP0-I , 1 ) 

35)  ERR0R=1  ! J . 0  * (SUM-1 . 0 ) 

Compute  percent  error  of  solution. 

36)  PRINT  9,T,ERR0R 

Print  current  time  T  and  percent  error. 

37)  9  F0RMAT ( 1 H  F15.4,F20.6) 

/  38)  11  RETURN 

Any  problem  statistics  such  as  the  expected  number  in 
the  system  would  be  computed  in  entry  point  NEWVAL  by  the 
user  as  required  by  his  particular  problem. 

WRAPUP  Entry  Point 

This  entry  point  is  called  from  the  Main  program  after 
the  problem  has  been  solved,  that  is,  when  time  T  reaches 
TMAX.  If  the  user  executes  a  RETURN,  the  Main  program  calls 
INIT,  thus  providing  the  user  the  opportunity  to  initialize 
a  new  problem.  A  CALL  EXIT  statement,  written  instead, 
terminates  the  run. 
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/  39)  ENTRY  WRAPUP 

40)  RETURN 

In  this  example,  the  RETURN  causes  a  call  to  INIT.  The 
READ  statement  is  executed  to  get  new  problem  parameters. 
If  no  more  data  is  in  the  reader,  a  branch  is  made  to 
statement  10  to  terminate  the  run. 

41)  10  CALL  EXIT 

/  42)  END 

The  FORTRAN  listing  of  this  program  and  sample  output 
are  included  next.  The  section  following  the  listings  will 
discuss  in  detail  the  problem  of  storage  of  the  probability 
vectors  in  the  P-array. 
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Dimensioning  the  Problem 


The  P-array  in  C0MM0N  storage  contains  the  solution 
vector  and  the  derivatives.  The  dimension  of  the  array, 
denoted  by  *  below,  is  constructed  to  correspond  to  the 
possibly  several  subscripts  of  the  user  problem  according  to 
normal  FORTRAN  conventions.  For  each  probability  there  is  a 
corresponding  derivative  with  the  same  subscript.  By 
appending  an  additional  subscript  on  the  P-array  with  the 
value  of  two,  two  identical  problem  areas  are  set  up.  The 
first,  that  is,  the  one  vrith  the  appended  subscript  value  of 
one,  is  for  the  solution  vector  probabilities.  The  second 
is  for  the  corresponding  derivatives.  The  compiler  assigns 
the  storage  in  ascending  order  so  that  the  first  half  is  the 
lower  part  and  the  second  is  the  upper  part.  Thus  P(*,2)  in 
the  C0MM0N  statement  sets  up  two  sequential  arrays  each  of 
dimension  *,  where  *  here  stands  for  the  dimension  of  the 
user  problem  as  discussed  below  in  more  detail.  The  user 
must  assign  NEQ  a  value  equal  to  the  product  of  the  sizes  of 
each  component  of  *  above.  The  main  program  is  then  able  to 
locate  the  derivative  locations  by  simply  adding  NEQ  to  the 
value  of  the  subscript  of  the  solution  vector  expression 
during  the  computations. 

The  Runge-Kutta  subroutine  requires  working  storage  for 
intermediate  results.  This  depends  on  the  size  of  the  user 
problem  and  thus  must  be  assigned  by  the  user  in  the  C0MM0N 
statement  as  the  W-array.  The  W-arrav  will  always  be 


dimensioned  as  exactly  one-and-one-half  times  the  size  of 
the  P  {*,2) -array.  An  example  will  oe  given  to  help  clarify 
this  point. 

Suppose  a  problem  has  equations 
Pijk(t)  =  f(t,Pijk(t)) 

where  i,j  and  k  are  indicies  used  to  enumerate  the  state  of  the 
systen ,  Suppose  further  for  this  nroblen  that 

1  <  i  £21 
i  1  j  <11 
1  <k  <6  . 

In  the  notation  used  above,  the  problem  dimension  *  is 
(21,11,6).  The  P-array,  however,  must  be  dimensioned  for 
this  problem  with  an  appended  subscript  of  value  two. 
Therefore,  in  order  to  make  space  available  for 
corresponding  derivative  values,  write  P  (21 , 1 1 ,6 ,2) .  Now 
the  W-array  must  be  one-and-one-half  the  size  of  the  P- 
array.  Write  W(4158),  which  is  the  product 
(21)  (11) (6) (2) (1.5)  . 

Occasionally,  a  problem  may  be  written  which  has  more 
subscripts  than  the  compiler  will  allow,  (For  the  IBM  360 
FORTRAN-IV  compiler  the  maximum  is  seven.)  Suppose  the 
problem  appears  as 
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(t)  . 


There  are  p  subscripts.  Also  suppose  these  run  from  one  to 
. ..,mp,  respectively.  The  p  subscripts  must  be 
replaced  by  a  lower  number  of  subscripts  for  the  compiler. 
Say  one  is  used,  denoted  by  j.  Given  a  set  of  p  indicies  in 
the  program,  the  following  algorithm  should  be  used  to 
compute  j  which  will  be  then  used  as  the  subscript  of  the  P- 
array . 

Initially  set  j=i  -1.  Then  compute 


j-«-jm  +i  -1  for  k=1  to  p-1, 
p-k  p-k 

and  then  compute 


j«-j  +  1. 


To  illustrate,  suppose  the  four  subscripts  of  the 

former  example  could  not  be  dimensioned  as  P{21, 11,6,2)  as 

was  done  because  of  compiler  restrictions.  To  use 

equivalent  single-subscripting  first  dimension  the  P-arrav 

for  the  product  of  the  subscripts,  P(2772).  Hero  m  =21,  m 

1  2 

=11,  m  =6,  m  =2  and  p=4.  To  locate  the  correct  subscript 
for  the  P-array  for  some  particular  choice  of  subscripts, 
say  (i  ,i  ,i  ,i  )=(2,4,5,1),  use  the  above  algorithm. 

12  3  4 

j=i  -1*1 -1=0,  then  for  k=1  to  3; 

H 


Set 
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j«-jm  +i  “1  =  (0)  (6)  +5-1=4 , 

t  i 


j+jm  +i  -1= (4) (11) +4-1*47, 
2  2 


j«-jm  +i  -1-(47) (2 1 ) +2-1=988, 


and  finally 


j«*j  +  1  =  988+1  =  989. 

Thus  P(989)  corresponds  to  P(2,4,5,1). 

Program  Output 

The  programming  package  treats  the  P-array  as  if  it 
were  singly-subscripted.  During  a  problem  solution,  the  P- 

array  will  be  printed  after  every  NPRNT  integration  steps  on 

» 

FORTRAN  unit  8.  The  output  will  consist  of  the  problem  NAME 
at  the  top  of  the  page  with  the  left  column  labeled  for  time 
T.  After  each  value  of  T  printed,  the  P-array  will  be 
printed  immediately  to  the  right.  Values  of  the  array  which 
are  less  than  5*10  will  not  be  printed  if  they  follow 
larger  values.  If  the  single-subscript  output  form  of  the 
P-array  is  not  desired  by  the  user,  define  FORTRAN  data  set 
8  as  DUMMY  in  the  operand  field  of  the  DD  control  card  for 
FT08F001  as  shown  below  in  the  Job  Control  Cards  section  and 
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Job  Control  Cards 

The  job  control  language  (JCL)  statements  for  running 
on  the  IBM  360/65  computer  at  the  University  of  Iowa  are  now 
presented. 

//QUEUEING  J0B  (XXXXXXXX) , ' STU  QUEUES* 

//  EXEC  PDQF0R , PARM=  *  LIST/PRINT  r MAP ' 

//S0URCE  DD  * 

(Place  source  deck  of  programming  package  and  subroutine 
INIT  here.) 

/* 

//DATA  DD  * 

(Place  any  data  to  be  read  by  subroutine  INIT  from 
unit  5  here.) 

/* 

//FT08F001  DD  S YS0UT= A , DCB=RECFM=UA 

(To  get  P-array  output  from  Main  program. ) 

//FT08F001  DD  DUMMY 

(To  ignore  output  from  Main  program.) 

//SYSABEND  DD  SYS0UT=A 

(Include  if  a  dump  is  desired  if  abnormal  termination 
such  as  OCU  or  0C5  occurs.) 


